High-accuracy virtual testing of air conditioner’s digital twin focusing on key material’s deformation and fracture behavior prediction

The concept of digital twin has been introduced for some time, yet one fundamental element of digital twin, digital material, has not been thoroughly studied. To interact with the physical product, the digital twin should always truthfully reflect the responses under various stimuli. In this paper, the deformation and fracture behavior of high impact polystyrene (HIPS) under the influencing factors of strain rate and stress triaxiality are studied to construct the material’s digital model. A digital twin of air conditioner product is further built and tested under virtual drop test. Comparing to experimental results, the acceleration curve, crazing induced whitening and the fracture events can all be captured by the digital twin. Our work demonstrates the importance of material characterization as an essential step to construct an accurate digital twin and shows a promising future of digital twin in virtual testing to replace traditional “trial and error” experiments.


Scientific Reports
| (2022) 12:12432 | https://doi.org/10.1038/s41598-022-16511-w www.nature.com/scientificreports/ HIPS is chosen as the model material to illustrate the importance of digital material for the accurate characterization of the digital twin. In domestic appliance industry, HIPS is often employed as outer casing acting as both high gloss appearance and structural support. After standard tests such as drop test, impact test and random vibration test, the product should not only retain its functions but also be free of whitening or visual cracks. This necessitates the building of the constitutive models of the involved materials such that a DT of the physical product could be constructed and tested virtually under various testing conditions. Due to its dual roles, the accuracy of the constitutive and failure model of HIPS becomes key to predict the mechanical behavior of the DT of domestic air conditioner at product development, production and operation stages. Through recent decades, the benefits of rubber content on improving the mechanical performance of HIPS has been thoroughly studied [23][24][25] . A 5-10% weight ratio of rubber content could typically yield a 20-30% volume content of rubber through phase inversion during polymerization 24,25 due to the formation of a special cellular or occluded structure, composed of many polystyrene inclusions embedded inside rubber 31 . Distinct from pure polystyrene, inside which concentrated crazes quickly connect to form critical crack and fracture, the rubber content in HIPS acts as both initiation sites and obstacles to craze; the latter prolongs the process for crazes to develop into critical crack length [32][33][34] . From microscopic view, crazes are commonly formed at equatorial zones of the rubber and propagate in direction perpendicular to tension. When the volume content of rubber exceeds about 15%, the close distance between rubber particles induces stress superposition, which allows thicker and longer craze to absorb more energy. The rubbers could further stop the craze from forming into long crack through crack tip blunting and force the crazes to form into distributed fashion enabling an enhancement of the ductility and impact resistance 25 . This phenomenon is widely observed under standard tensile tests of HIPS [35][36][37][38] . From macroscopic view, it has been shown that the threshold stress for crazing in HIPS is lower than the stress required for shear yielding 39 in tension and the toughness under tension mode (mode I) is typically lower in comparison to that under shear mode (mode II) 40,41 . Thus, many researchers believe under most load-bearing scenarios the final fracture of HIPS should be induced by tension and accompanied with features of crazing 42,43 . However, there still lacks some evidence on how stress state other than tension would affect the fracture modes. Besides, the fracture pattern could also be altered based on different craze-crack interactions 24,25,[44][45][46][47] . Higher loading rate would typically cause the craze opening rate to accelerate 48 and altering the fracture pattern from features like parabola shapes at low crack extension speed to broken patch patterns at high crack extension speed; this phenomenon is found for both polystyrene 24 and polypropylene 45 ; whether similar phenomena could be found for HIPS remains unknown. Therefore, how factors like strain rate, stress state and fracture pattern could affect the deformation and fracture behavior of HIPS is still not fully explored and have become obstacles to construct detailed constitutive and failure model of HIPS, which are crucial for accurate virtual testing of the DT.
In this paper, the mechanical behavior of HIPS under the strain rate range of 0.001-100 s −1 and various designed stress triaxiality will be examined. Fractographic analysis of fracture surfaces of chosen stress state is conducted to further understand the damage-fracture mechanism of HIPS. A phenomenological constitutive and failure model of HIPS is constructed and incorporated in a digital twin of an air conditioner product with HIPS as outer casing and base plate. The product and its digital twin are further tested under drop test experimentally and virtually, respectively and compared to demonstrate the importance of digital materials in high accuracy prediction for digital twin.

Methods
Structure of the digital twin. In order to realize the possible deformation and fracture behavior of air conditioner product at various mechanical testing conditions, it is necessary to establish a digital twin system to truthfully reflect its responses. The main approach employed here is "bottom-up" as shown in Fig. 1: i.e. start from digital material to digital twin and physical product. At digital material level, the deformation and fracture behavior of HIPS is studied. Through characterization and analysis, the digital model of material could be constructed. At digital twin level, each part/component is dawn by computer aided design (CAD), converted, combined with material digital models and assembled into product through computer aided engineering (CAE). The digital twin is further verified under drop test in a consecutive manner. Quantitative comparisons are made based on experimentally collected data from physical product and data predicted by its DT. Qualitative comparisons are made based on crazing induced whitening, fracture path and fracture pattern observed after physical product test and the corresponding phenomena predicted by the DT directly or indirectly. At physical product level, real time data is captured by the product such that real time predication on the product deformation could be made to help identify the product quality and increase efficiency of the decision making process.
Material and testing. Plastic sheets of HIPS are cut directly from the injection-molded outer casing of product and then machined through computer numerical control (CNC) technique to specific geometry of test specimen. Figure 2 presents the 7 specimen designs; for static and dynamic tensile test, the specimens are designed in accordance with ISO 527 and ISO 8256, respectively; for other tests the specimen designs are customized to include different stress states. The detailed types of experiments are listed in Table 1. The ideal stress triaxiality covers the range of − 0.33 to 0.66. The static tensile test is conducted using universal tensile machine (UTM) Instron 5982 with strain rate of 0.001 s −1 . The dynamic tensile test is conducted using ZwickRoell HTM 5020 high-speed testing machines with strain rate ranging from 0.1 to 100 s −1 . All the strain data is obtained through digital image correlation (DIC) technique. At least three repeated experiments of each type are conducted to ensure the accuracy of the experiments. The fracture path is observed by Keyence optical microscope (OM) and the fracture surface is observed using scanning electron microscope (SEM) Tescan VEGA3. The acceleration of the product during drop test is measured using Brüel & Kjaer accelerometer sensors and analyzed using its dynamic signal analyzer. Our method focuses on building the digital twin of air conditioner that can predict the mechanical responses of physical product under mechanical stimuli. The prototype DT is tested under drop test without cushion protection to mimic an extreme condition while being able to be tested under our testing facilities and collect data to be compared. In addition, two consecutive tests of the DT allow us to consider the effect of loading history on accumulative damage of the product, which typical simulation would not consider.  To construct the materials' constitutive models, the hardening curve under different strain rate conditions are extracted from true-stress strain curve, fitted to proper formula with accuracy over 90%, extrapolated to 150% strain and incorporated in piecewise linear plasticity material model.
The failure model employed is generalized incremental stress-state dependent damage (GISSMO) model, which consider stress softening and accumulated damage to correlate failure events 49 . The key input is equivalent plastic strain to failure curve constructed based on the iterative algorithm of CrachFEM, a commercial software. The primary goal is to minimize the residual error between the calculated results of the finite element model and the force displacement curves from actual tests 50 . The digital twin of the air conditioner is constructed with all sub-systems and associated components as shown in Fig. 3a-c. The element size is optimally chosen as 3 × 3 mm after mesh convergence study. The overall model consists of 2 accelerometer element, 302,146 shell elements and 492,457 solid elements with a total number of 794,605 elements. Two consecutive drop tests at product level are designed. The first drop test corresponds to the vertical surface drop test, which aims to collect the acceleration data as the main parameter to verify the accuracy of the DT. Two accelerometer sensors are placed on the top and bottom of the physical product as shown in Fig. 3c-g. The second drop test is of edge drop type with a predefined inclination angle in order to create apparent crack and stress whitening such that the failure model can be verified. It should be noted the area of impact is meshed with smaller element i.e. 1 × 1 mm to better reflect the detailed crack path. A strain mapping technique is also employed to account for the plastic strain history from 1st surface drop test into 2nd edge drop test.

Results
Static and dynamic tensile test results. The static and dynamic responses are fundamental to understand the mechanical behavior of HIPS. The stress strain curves of HIPS are shown in Fig. 4. The strain hardening rate curve is also extracted and analyzed. 3 deformation stages can be identified: stage I as the elastic region corresponding to strain level of 0-0.02, stage II as strain softening region with 30% reduction in stress corresponding to strain level of 0.02-0.1 and stage III as strain hardening region corresponding to strain level of 0.1 up to the strain at final fracture.
The stress strain curves of HIPS under strain rate ranging from 0.001 to 100 s −1 are shown in Fig. 4b. It can be observed as strain rate increases, the yield strength increases significantly whereas the ductility drops considerably due to the loss of stage III. Constitutive, failure model and validation. In order to construct accurate constitutive and failure model, some key parameters required for FE analysis need to be identified. The elastic modulus is estimated by linear fit to the experimental data at all strain rates as shown in Fig. 5a and calculated as 2041 ± 113 MPa. It can be seen that the strain rate has almost no effect on the elastic modulus. It is found the hardening model can be constructed by semi-experimental and semi-fitting method 51 , i.e. the experimental data up to stage II is kept and a quadratic polynomial fit is added and extrapolated to 100% strain as shown in Fig. 5b. The plastic strain failure curve for input of GISSMO failure model is extracted from CrachFEM and fitted with quantic polynomial as shown in Fig. 5c.
By putting the necessary parameters into the constitutive and failure model, a complete digital model of HIPS is constructed. To validate the model, the load-displacement curves obtained through 6 designed tests are compared to those calculated from simulations, respectively as shown in Fig. 6. It can be seen the results obtained under simulation are consistent with the experimental results at various tested stress state conditions. Fracture surface observations. To understand the damage-fracture mechanism of HIPS, the fracture surfaces from chosen tests are observed. The fracture surface of HIPS under shear test (Fig. 7) shows two completely distinct regions: region A has very fine features whereas region B has large fibrous sheets, which extend fully to the thickness of the specimen. Crack propagation lines are clearly visible as indicated by the white arrows shown inside region A in Fig. 7b. Inside region A, slanted torn features and numerous tiny voids are clearly seen during crack extension (Fig. 7c,d). Inside region B, the large fibrous sheets are resultant from the final instantaneous fracture (Fig. 7a,c), which was confirmed by the video recording.
Under both tensile test and notched tensile test, typical ductile fracture morphology is observed for HIPS. At low magnification, Fig. 8a,b indicates both samples illustrate crack interactions, whereas a ridge is formed when    Virtual and experimental results of digital twin and physical product. To further verify the digital model of HIPS, the DT of an air conditioner is constructed, tested virtually and compared with physical product's experimental results. Figure 9a illustrates the acceleration curves obtained from both virtual and experimental surface drop tests. For the first peak of acceleration, the difference between the digital twin and the physical product for both top and bottom accelerometer sensors are quite small. Interestingly, the consistency of second peak for the bottom is better than the top, which might be due to the position of the sensor, where a stiffer component is placed on the bottom (Fig. 3g) than the top (Fig. 3e). The energy map of the digital twin, Fig. 9b indicates that overall, the total energy is conserved and the hourglass energy is close to 0. The kinetic energy from free fall drop is converted to internal energy of the system of the digital twin. From the perspective of acceleration prediction and energy conservation, the DT is verified. The second drop test is of edge drop type with a pre-defined inclination angle of 30° and height of 1 m between the impact edge and the ground. Apparent crack and stress whitening can be observed after test. Figure 10 illustrates the overall crack profile after actual edge drop test and that predicted by the DT. The major crack extension and crack initiation sites are well captured by the DT as shown in Fig. 10a,b. To better show the crack morphology, the part containing the major crack is cut and viewed under OM (Fig. 10d,e). Figure 10c shows the matching part containing this crack from the DT; it can be seen the detailed crack deviation points and path developments of the DT are all in agreement with those that appears in physical product (Fig. 10d,e). Furthermore, strands of fibrous sheets, which extend fully to the thickness is identified as type 1 fracture as shown in Fig. 10d.
To further understand the fracture mechanism, the crack initiation site corresponding to the highest stress concentration is examined under SEM. Very different fracture features are observed as shown in Fig. 11. At higher magnification, it can be seen type 2 fracture surface is composed of typical mackerel or patch morphology (Fig. 11b), type 3 fracture surface contains highly plastically stretched fibrils with sharp tip and strong orientations (Fig. 11c) and type 4 fracture surface contains rugged surface with majority of half or unfinished circular steps (Fig. 11d).

Discussions
Deformation behavior of HIPS. Under static tensile test, 3 deformation stages can be defined for HIPS, similar to polymers like PMMA, PC-ABS alloy 51 and polycarbonate (PC) 52 . However, the mechanisms behind for those materials are quite different, this can be seen by the different necking developments and strain hardening rate curves. For HIPS and PMMA, crazing is the dominant deformation mechanism causing the initial necking to develop uniformly and total volume to increase before final fracture (Fig. 12a,b); since the strain extension is balanced by crazing only the strain hardening rate is quite low (Fig. 4a). For PC, ABS and PC-ABS alloy, the deformation is controlled by the competition between crazing and shear yielding such that apparent necking or strain localization (Fig. 12c) can be observed before fracture; shear yielding, which conserves the volume similar to strain hardening stage of metals, also causes the higher strain hardening rate, higher strength and lower ductility than those of HIPS 51 . This finding suggests that despite the price advantage of HIPS comparing to ABS, a reconfiguration prediction of the DT before material replacement is necessary due to their quite different mechanical properties.
With increasing strain rate, Fig. 4b shows the strength increases yet the ductility deceases with the disappearance of craze development stage, stage III. Since the resistance to polymer chain alignment causes strain hardening in the material at high strain rate 53,54 , this process can also be considered an energy-activated and rate-dependent process, an equation based on Eyring's equation is expressed as 55 : where σ y = yield strength, R = ideal gas constant, T = temperature, v = activation volume, ė = strain rate, ė 0 = constant pre-exponential factor, ΔH = height of potential energy barrier and k = Boltzmann's constant. This equation implies that there exists a linear relationship between σ y and ln(ė). Such linear relationship is also found for the present HIPS (Fig. 12d), in agreement with findings on polymers like poly(viny1 chloride) PVC 56 . (e) (f) (c) (d)  www.nature.com/scientificreports/ It is speculated that the decrease in ductility at high strain rate might be related to 2 reasons: the increase in the density of micro voids or cavitation due to the less time for chain disentanglement 54 and the possibility that rubber could not relax to build up enough stress to withstand the load with matrix together 55 . The increment of strength and altered hardening behavior were both taken into account in the material model, which helps in identifying crazing induced whitening as show in "Deformation induced visual whitening prediction".
Effect of stress triaxiality on fracture behavior of HIPS. The effective plastic strain to failure as a function of stress triaxiality (shown in Fig. 5c) is extracted by CrachFEM's internal algorithm, which considers both normal fracture and shear fracture 50 . The curve is drawn through the iterative calculation results aiming to minimize the residual error between the calculated results of the FE model and the actual test force displacement curves, it can be seen the curve is very similar to the normal fracture curve indicating that at stress triaxiality greater than 0, thin sheet HIPS with plane stress states tends to fail by normal fracture mode, i.e. coalescence of micro-voids instead of shear band development for shear fracture. However, it seems the initial stress triaxiality and crack propagation speed could still influence the fracture pattern, our observation on the fracture surfaces of HIPS under different testing conditions brings some new insight on this issue.
A fracture pattern map of HIPS resultant from the interaction of initial stress triaxiality, η and crack propagation speed, V c is proposed and shown in Fig. 13. At quadrant of low η (shear mode) and high crack propagation speed, type 1 fracture surface contains large fibrous sheets (Fig. 7c, region B). It is speculated this morphology is due to the later tensile rupture of the shear band that formed during the beginning of shear test. Very similar phenomenon was observed in Ref. 57 , which subject atactic polystyrene to compression first and then tension. At quadrant of high η (0.33 and 0.44) and high crack propagation speed, type 2 fracture surface is very smooth showing patch patterns, which is consistent with findings on the fracture surfaces observed under high crack extension rate 39 and low temperature 58 . At quadrant of low start η (shear mode) and low crack propagation speed (static condition), the fracture surface is of type 3, composed of very fine, highly stretched and fibrillated material (Fig. 7d) indicating that even at shear-dominate test the fracture is still controlled by craze-crack interaction at fine scale. It should be pointed out that this is due to the inability of actual test to maintain a pure shear condition www.nature.com/scientificreports/ and the stress triaxiality would change to be positive as test continues. Similar morphology can be seen in Ref. 25 . At quadrant of high η and low crack extension speed, the fracture surface is of type 4 showing half-circle pattern formations that are found after tensile and notched tensile test (Fig. 8), which might be resulted from interactions between growing circular crazing and a moving planar crack.
Validation of air conditioner digital twin. The digital material model of HIPS is further utilized and validated in virtual drop test of the DT of an air conditioner product. The digital twin has been verified first by the correlation of the acceleration curves between predicted data and that measured by accelerometer and the www.nature.com/scientificreports/ energy conservation from surface drop test as shown in "Virtual and experimental results of digital twin and physical product". The detailed deformation and fracture events predicted by the DT in comparison to those that occurred in physical product at second inclined drop test are shown in the following: Deformation induced visual whitening prediction. For HIPS, it has been proved that crazing is the main deformation mechanism under most loading scenarios, thus the plastic strain should be an important indicator of the crazing induced whitening. Based on this understanding, three crazing criteria are established to predict the crazing induced whitening. 0.025, 0.05 and 0.1 plastic strain are chosen as the crazing criteria separately and the respective plastic strain contours of the DT are drawn for the two parts that experience major impact. For the base plate, in comparison to the whitening area after the drop test (Fig. 14b), it is shown that under both criteria of 0.025 (Fig. 14c) and 0.05 (Fig. 14d) the predicted whitening areas (shown in red) are overestimated whereas only under the 0.1 criterion both the locations and sizes of the whitening zones are well captured, as indicated by the white arrows in Fig. 14b,e. Similar trend can also be seen for the observed area on the outer panel shown in Fig. 15. It is worth mentioning; the 0.1 plastic strain criterion is in accordance with the starting point of stage III (necking development) shown in Fig. 4a, correlating the accurate prediction of our DT to experimentally measurable material data.
Fracture prediction. As shown in "Virtual and experimental results of digital twin and physical product", sound correlations can be established on the crack origin and crack extension profiles between the results predicted by the DT and that observed on the fractured part from physical product under OM. The upper fractured part (Fig. 10c) predicted by the DT matches the real fractured part (Fig. 10d,e) at majority of the crack-developing sites. In addition, the stress triaxiality values of those elements corresponding to the crack initiation site (Fig. 16a) are drawn in Fig. 16b. It can be seen the value of stress triaxiality changes in 2 stages: in the first stage it decreases to negative, increases to positive and then becomes stable at high stress triaxiality value in the second stage. Such varying stress state suggests that mixed modes of fracture might exist, which are confirmed by the different fracture patterns observed at crack initiation site under SEM (Fig. 11).   Fig. 16c for comparison. It can be seen that at crack initiation site the change of stress triaxiality in the first stage is in correlation to that of individual shear test i.e. the stress triaxiality value under both curves first decreases to negative value and then increases to positive value. This explains the type 3 fracture surface observed under both cases, which contains highly plastically stretched fibrils. In the second stage, the stress triaxiality curve becomes flat at high stress triaxiality value, which suggests a combination of tensile and notched tensile. Indeed, type 2 patch patterns and type 4 unfinished circular step patterns are observed at and near crack origin, respectively (Fig. 11) corresponding to the fracture surfaces observed at individual tensile and notched tensile test (Fig. 8). It is worth mentioning away from crack origin, Fig. 10 shows a portion of the fracture surface along the crack path is composed of large fibrous sheets, which corresponds to type 1 fracture pattern. According to 57 , the fibrous sheets inside atactic polystyrene are caused by tension of the shear band that formed previously under compression. This is also in agreement with our findings on the final fracture surface of HIPS after shear test (Fig. 7 region B), which resulted from the shifted stress state from pure shear to tensile and high crack extension rate.
Our work indicates that accurate digital material models could enhance the DT of consumer products at all 3 phases. At design phase, the feasibility of materials to be employed in components of consumer products could be evaluated in various aspects. Local deformation, fracture events of parts and cushioning property (acceleration) of product could all be captured as crucial guidance for adjustments of product design parameters. At reconfiguration stage, different product structure, materials and connection types could all be evaluated. Especially, due to the different deformation behaviors of HIPS and ABS, it is suggested a reconfiguration predication of the DT is necessary before material replacement. A physical product of the air conditioner is further implemented into drop test. High level of accuracy of the DT is achieved both quantitatively and qualitatively. In particular, the technique of microscopic characterization of the fracture surfaces is employed as an innovative verification method for testing the accuracy of the DT. At operation stage, based on the input of the physical product, the www.nature.com/scientificreports/ DT would predict and give instructions to the physical product. In our case, the accelerometer readings from physical product are treated as input to the DT, the DT could predict the material behavior of all parts and make judgments based on associated standards. It should be pointed out that currently our research focus lies in the accurate prediction of the mechanical and fracture behavior of the DT only within the scenario of drop test, which limits a deeper interaction between our digital twin and product. In the next stage, the DT would be tested under other tests including inclined impact, random vibration and compression tests to verify and improve the accuracy of the DT. In addition, due   www.nature.com/scientificreports/ to currently limited computing speed for simulation, the mechanical responses of the DT under various testing conditions of each test should be recorded and related to physically measurable mechanical stimulus such that a deformation database of the DT could be constructed and retrieved. For the physical product, wireless sensors would be installed on the product package to reflect the mechanical stimuli in "real-time" such that its the DT could reflect the deformation and damage concurrently to judge if the product passes the associated standards.

Conclusion
In this paper, the deformation and fracture behavior of HIPS is investigated. The influences of strain rate and stress triaxiality are considered. The newly obtained knowledge is further incorporated in the material model and utilized to predict the detailed whitening and fracture phenomena of an air conditioner's digital twin through virtual testing. From material's perspective, it is found the deformation behavior of HIPS under tension consists of three stages: elastic, strain softening and strain hardening. The uniform necking development up to final fracture observed under DIC suggests that crazing is the dominant mechanism, which differs from ABS with local necking before fracture. With increasing strain rate, the ductility is greatly reduced and yield strength increases linearly with natural logarithm of strain rate. In addition, although the fracture surfaces are resultant of crazing and crack front, the stress triaxiality and crack propagation speed are identified as two controlling factors for final fracture patterns. A fracture map is constructed based on present findings.
From product's perspective, by incorporation materials' digital models into the digital twin of an air conditioner, its deformation and fracture behavior under drop test can be predicted. For deformation prediction, a crazing criteria of 10% plastic strain, identical to the beginning of crazing development stage in tensile test, is identified and used to accurately predict both locations and degrees of the whitening for outer casing and base plate of the DT in agreement with those of physical product. For fracture prediction, a fracture model based on plastic failure strain and stress triaxiality curve is employed in the DT. The fracture path captured by the DT is also observed under OM at the same location of physical product. In addition, by extracting the stress triaxiality history of the fracture initiation site of the DT and relating to the fracture pattern map, mixed fracture patterns are predicted. This is verified by the fracture surface observation of the corresponding site on physical product under SEM, which consists of finely ruptured craze, patch pattern and half-circle pattern resultant from the change of stress state. Our work demonstrates the importance of material characterization for the high accuracy of the digital twin of consumer products and shows a promising future of the digital twin to replace traditional "trial and error" experiments, monitor "real-time" product life cycle and achieve predictive maintenance.

Data availability
The datasets generated during and/or analysed during the current study are not publicly available due to ongoing research but are available from the corresponding author on reasonable request.